﻿/*
 * Created by SharpDevelop.
 * User: Administrator
 * Date: 2016/11/11
 * Time: 16:51
 * 
 * To change this template use Tools | Options | Coding | Edit Standard Headers.
 */
using System;
using MapProj;
namespace GeoFly
{
	/// <summary>
	/// Description of Sampling.
	/// </summary>
	public static class Sampling
	{
		private static int Search(double value, double[] array)
		{
			int len1 = array.Length;
			if (value <= array[0])
				return 0;
			if (value >= array[len1 - 1])
				return len1 - 1;
			for (int i = 0; i < len1; i++) {
				if (value < array[i])
					return i - 1;
			}
			throw new Exception("Invalid result");
		}
		private static int SearchReverse(double value, double[] array)
		{
			int len1 = array.Length;
			if (value >= array[0])
				return 0;
			if (value <= array[len1 - 1])
				return len1 - 1;
			for (int i = 0; i < len1; i++) {
				if (value > array[i])
					return i - 1;
			}
			throw new Exception("Invalid result");
		}
		private static double LinearInter(double AValue, double BValue, double ACoor, double BCoor, double ECoor)
		{
			if (BCoor - ACoor == 0)
				return (AValue + BValue) / 2;
			return (AValue * (BCoor - ECoor) + BValue * (ECoor - ACoor)) / (BCoor - ACoor);
		}
		/// <summary>
		/// Interpolate data using bilinear functions
		///        newlats,newlons is the destination coordinates
		///			oldlats,oldlons is the original coordinates of data
		///    example(for the domain of China): 
		///        ReGrid(data,numpy.arange(0,61),numpy.arange(70,141),lat,lon,ReverseLat=True)
		/// </summary>
		/// <param name="newlat"></param>
		/// <param name="newlon"></param>
		/// <returns></returns>
		public static double[,] ReGrid_NotEven(double[,] data, double[] latsSource, double[] lonsSource, double[] newlats, double[] newlons, bool ReverseLat)
		{
			int nlat = newlats.Length;
			int nlon = newlons.Length;
           
			double[,] result = new double[nlat, nlon];

			int[] lonIndice = new int[nlon];

			for (int i = 0; i < nlon; i++) {
				double lon = newlons[i];
				lonIndice[i] = Search(lon, lonsSource);
			}
			for (int ilat = 0; ilat < nlat; ilat++) {
				double lat = newlats[ilat];
				int latIndex = -1;
				if (ReverseLat == false) {
					latIndex = Search(lat, latsSource);
				} else {
					latIndex = SearchReverse(lat, latsSource);
				}
				for (int ilon = 0; ilon < nlon; ilon++) {
					//Calculate bilinear interpolation using nearest four grid cells
					double lon = newlons[ilon];
					//lonIndex=Search(lon,lons1);
					int lonIndex = lonIndice[ilon];

					double A = data[latIndex, lonIndex];
					double B = data[latIndex, lonIndex + 1];
					double C = data[latIndex + 1, lonIndex];
					double D = data[latIndex + 1, lonIndex + 1];
                    
                  
                    	

					double lonA = lonsSource[lonIndex];
					double lonB = lonsSource[lonIndex + 1];
					double E = LinearInter(A, B, lonA, lonB, lon);
					double F = LinearInter(C, D, lonA, lonB, lon);

					double latA = latsSource[latIndex];
					double latB = latsSource[latIndex + 1];
					result[ilat, ilon] = LinearInter(E, F, latA, latB, lat);
				}
			}
			return result;
		}
		public static double[,] ReGridBilinear(double[,] data, double[] latsSource, double[] lonsSource, double[] newlats, double[] newlons)
		{
			int nlat = newlats.Length;
			int nlon = newlons.Length;
           
			double[,] result = new double[nlat, nlon];

			int[] lonIndice = new int[nlon];
			double dlat = -(latsSource[1] - latsSource[0]);
			double dlon = lonsSource[1] - lonsSource[0];
			if(dlat<0)
			{
				Array.Reverse(latsSource);
				dlat*=-1;
			}
			for (int i = 0; i < nlon; i++) {
				double lon = newlons[i];
                
				lonIndice[i] = (int)((lon - lonsSource[0]) / dlon);// Search(lon, lonsSource);
			}
			for (int ilat = 0; ilat < nlat; ilat++) {
				double lat = newlats[ilat];
				int latIndex = -(int)((lat - latsSource[0]) / dlat);// Search(lat, latsSource);
				//latIndex=latsSource.Length-1-latIndex;
				for (int ilon = 0; ilon < nlon; ilon++) {
					//Calculate bilinear interpolation using nearest four grid cells
					double lon = newlons[ilon];
					//lonIndex=Search(lon,lons1);
					int lonIndex = lonIndice[ilon];

					double A = data[latIndex, lonIndex];
					double B = A;
					if (lonIndex + 1 < lonsSource.Length)
						B = data[latIndex, lonIndex + 1];
					double C = A;
					double D = A;
					if (latIndex+1 < latsSource.Length) {
						C = data[latIndex+1, lonIndex];
						if (lonIndex + 1 < lonsSource.Length)
							D = data[latIndex+1, lonIndex + 1];
					}
					
					double[] abcd=new double[]{A,B,C,D};
					bool flag=false;
					for(int ind=0;ind<4;ind++)
					{
						if(abcd[ind]<-9990)
						{
							result[nlat- ilat-1,ilon]= Math.Max(A,Math.Max(B,Math.Max(C,D)));
							flag=true;
						}
					}
					if(flag)
						continue;
					double lonA = lonsSource[lonIndex];
					double lonB = lonA;
					if (lonIndex + 1 < lonsSource.Length)
						lonB = lonsSource[lonIndex + 1];
					double E = LinearInter(A, B, lonA, lonB, lon);
					double F = LinearInter(C, D, lonA, lonB, lon);

					double latA = latsSource[latIndex];
					double latB = latA;
					if (latIndex + 1 < latsSource.Length)
						latB = latsSource[latIndex+1];
					double vvv= LinearInter(E, F, latA, latB, lat);
					result[nlat-1- ilat, ilon] =vvv;					
				}
			}
			return result;
		}
		
		private static double Interp(double curLat,double curLon,double[,] data,double[] latsSource, double[] lonsSource)
		{
	        double dlat=-(latsSource[1] - latsSource[0]);
	        double dlon = lonsSource[1] - lonsSource[0];
			int latIndex= -(int)((curLat - latsSource[0]) / dlat);
			int lonIndex= (int)((curLon - lonsSource[0]) / dlon);
			if(dlat<0)
			{
				
				dlat*=-1;
			}
			if(lonIndex+1>=lonsSource.Length)
			{
				return -9999;
			}
			double A=data[latIndex  ,lonIndex];
			double B=data[latIndex,lonIndex+1];
			double   C=data[latIndex+1,lonIndex];
			double D=data[latIndex+1,lonIndex+1];
			double lonA=lonsSource[lonIndex];
			double lonB=lonsSource[lonIndex+1];
			double	E=LinearInter(A,B,lonA,lonB,curLon);
			double  F=LinearInter(C,D,lonA,lonB,curLon);
			double  latA=latsSource[latIndex];
			double  latB=latsSource[latIndex+1];
			if(A<-9990 || B<-9990 || C<-9990 || D<-9990)
			return Math.Max(A,Math.Max(B,Math.Max(C,D)));
			return LinearInter(E,F,latA,latB,curLat);
		}
		
		public static double[,] ReGridLambert(double[,] data, double[] latsSource, double[] lonsSource, double startLat, double startLon,double cellSize,int rowCount,int colCount
		                                     ,double refLat,double refLon,double stdLat1,double stdLat2)
		{
			Array.Reverse(latsSource);
			double[,] result = new double[rowCount, colCount];
			MapProj.GaussKruger proj=new GaussKruger();
			double[] values=proj.ToXY_Lambert(startLon,startLat,refLon,refLat,stdLat1,stdLat2);
			double lux=values[0];
			double luy=values[1];
			double X0=0;
			double Y0=0;
			for(int row=0;row<rowCount;row++)
			{
				for(int col=0;col<colCount;col++)
				{
					X0=lux+cellSize*col;
					Y0=luy-cellSize*row;
					double[] vv= proj.ToLonLa_Lambert(X0,Y0,refLon,refLat,stdLat1,stdLat2);
					double lat= vv[1];
					double lon=vv[0];
					result[row, col]=Interp(lat,lon,data,latsSource,lonsSource);
				}
			}
			return result;
		}
		private static double AbsDist(double x0,double y0,double x1,double y1)
		{
			return Math.Abs(x0-x1)+Math.Abs(y0-y1);
		}
		public static double[,] ReGridNearest(double[,] data, double[] latsSource, double[] lonsSource, double[] newlats, double[] newlons)
		{
			int nlat = newlats.Length;
			int nlon = newlons.Length;
           
			double[,] result = new double[nlat, nlon];

			int[] lonIndice = new int[nlon];
			double dlat = -(latsSource[1] - latsSource[0]);
			double dlon = lonsSource[1] - lonsSource[0];
			if(dlat<0)
			{
				Array.Reverse(latsSource);
				dlat*=-1;
			}
			for (int i = 0; i < nlon; i++) {
				double lon = newlons[i];
                
				lonIndice[i] = (int)((lon - lonsSource[0]) / dlon);// Search(lon, lonsSource);
			}
			for (int ilat = 0; ilat < nlat; ilat++) {
				double lat = newlats[ilat];
				int latIndex = -(int)((lat - latsSource[0]) / dlat);// Search(lat, latsSource);
				//latIndex=latsSource.Length-1-latIndex;
				for (int ilon = 0; ilon < nlon; ilon++) {
					//Calculate bilinear interpolation using nearest four grid cells
					double lon = newlons[ilon];
					//lonIndex=Search(lon,lons1);
					int lonIndex = lonIndice[ilon];

					double A = data[latIndex, lonIndex];
					double ADist=AbsDist(lat,lon,latsSource[latIndex],lonsSource[lonIndex]);
					double B = A;
					double BDist=ADist;
					if (lonIndex + 1 < lonsSource.Length)
					{
						B = data[latIndex, lonIndex + 1];
						BDist=AbsDist(lat,lon,latsSource[latIndex],lonsSource[lonIndex+1]);
					}
					double C = A;
				   double	CDist=ADist;
					double D = A;
					double DDist=ADist;
					if (latIndex+1 < latsSource.Length) {
						C = data[latIndex+1, lonIndex];
						CDist= AbsDist(lat,lon,latsSource[latIndex+1],lonsSource[lonIndex]);
						if (lonIndex + 1 < lonsSource.Length)
						{
							D = data[latIndex+1, lonIndex + 1];
							DDist=AbsDist(lat,lon,latsSource[latIndex+1],lonsSource[lonIndex+1]);
						}
					}

					double[] aaa=new double[]{A,B,C,D};
					double[] bbb=new double[]{ADist,BDist,CDist,DDist};
					double min=1e20;
					double minValue=0;
					for(int i=0;i<4;i++)
					{
						if(min>bbb[i])
						{
							min=bbb[i];
							minValue=aaa[i];
						}
					}
					result[nlat-1- ilat, ilon] =minValue;				
					
				}
			}
			return result;
		}
		
	
	}
}
